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We present a numerical and partially analytical study of classical particles obeying a Langevin 
equation that describes diffusion on a surface modeled by a two dimensional potential. The potential 
may be either periodic or random. Depending on the potential and the damping, we observe 
superdiffusion, large-step diffusion, diffusion, and subdiffusion. Superdiffusive behavior is associated 
with low damping and is in most cases transient, albeit often long. Subdiffusive behavior is associated 
with highly damped particles in random potentials. In some cases subdiffusive behavior persists over 
our entire simulation and may be characterized as metastable. In any case, we stress that this rich 
variety of behaviors emerges naturally from an ordinary Langevin equation for a system described 
by ordinary canonical Maxwell-Boltzmann statistics. 
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I. INTRODUCTION 



Diffusion processes of atoms, molecules, and clusters 
of molecules on surfaces have been subjects of research 
for many decades due to their intrinsic interest and their 
technological importance. Some examples of modern ap- 
plications include self-assembled molecular filmgrowth, 
catalysis, and surface-bound nanostructures 0,01 ■ Also, 
many techniques that are used to characterize the growth 
of surfaces are based on the diffusion and subsequent ad- 
sorption of particles or molecules. These applications 
have led to a recent resurgence of interest in such pro- 
cesses, but now involving the motion of small and large 
organic molecules 0, 0] and of adsorbed metal clusters 
composed of tens and even hundreds of atoms @, @ . 

Recent research activity has been both experimental 
and theoretical QEHE&S0ES[Ifl[HEi, 
the principal questions in these studies being the de- 
termination of the jump lengths of large molecules or 
clusters on surfaces, and the diffusion asymmetry along 
some surface crystallographic directions or along partic- 
ular paths when there are many obstacles such as steps 
and/or impurities to negotiate. One of the more excit- 
ing experiments has involved the direct observation of 
the surface mobility of fairly large organic molecules us- 
ing a newly devised method of independently measur- 



ing the mean square displacement and the hopping rate 
of these molecules, and then taking advantage of a well 
known relation between them to extract the root mean 
square jump lengths j4|. These experiments seem to in- 
dicate that, as with single atoms in some cases [l3l Il4j . 
long jumps spanning multiple lattice spacings may play 
a dominant role in the diffusion of these molecules. An- 
other "experiment" , this one numerical (molecular dy- 
namics), has led to the prediction that clusters of hun- 
dreds of Au atoms will exhibit Levy-typepower-law flight 
length and sticking time distributions 6]. On the other 
hand, the theoretical literature on the subject of surface 
diffusion tends to be rather complex and, most impor- 
tantly, tends to invoke Levy walks or nights [l5| as a 
model input. While these approaches can provide in- 
sights on the effects of other elements of the model on 
the dynamics of particles on surfaces, they provide little 
insight on how the motion involving long jumps comes 
about in the first place. These simulations, which aim 
to reproduce realistic systems, do not stress the generic, 
minimal ingredients of such behavior. 

While a microscopically accurate analysis of surface 
diffusion requires extensive calculations (e.g., ab initio, 
or molecular dynamics), even the most powerful cur- 
rently available computers can not carry such calcu- 
lations to any where near experimentally relevant time 
scales [14l Il6| . Moreover, current experimental probes 
of the topography of surfaces, scanning tunneling mi- 
croscopy and atomic force microscopy, are usually carried 
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out at relatively high temperatures, which leads to addi- 
tional difficulties for first-principles calculations. There- 
fore, simpler approaches are essential and valuable [l7j . 
The scenario developed in this paper is of the simple 
Langevin type. It models classical particles moving in a 
two-dimensional potential, periodic or random, under the 
action of (5-correlated Gaussian thermal fluctuations and 
the associated linear dissipation, the important control 
parameter of the model being the friction coefficient. In 
our search for different behaviors as a function of the con- 
trol parameter, we must be mindful of the fact that exper- 
imental observations are not necessarily strictly asymp- 
totic, that is, that a particular behavior may be transient. 
Transient behaviors may of course persist for very long 
times, but to find them theoretically it is not enough to 
carry out an asymptotic analysis. We find that in spite 
of the simplicity of the model, it is able to reproduce 
the entire range of experimentally and computationally 
observed phenomenologies, ranging from superdiffusion 
through large-step diffusion all the way to subdiffusion. 
In particular, we show that Levy-like statistics appear 
quite naturally within the usual Langevin scheme for un- 
derdamped motion in a periodic potential. It is impor- 
tant to stress that in this case the Levy-behavior emerges 
as an intermediate asymptotic for an ensemble of parti- 
cles with a Maxwell-Boltzmann distribution of velocities. 

Our paper is organized as follows. In Sec.[n]we present 
the Langevin model and briefly list the quantities to be 
considered in subsequent sections to characterize the mo- 
tion of an ensemble of particles on a surface. In Sec. IIIII 
numerical results are presented for the periodic surface 
potential, as are analytic results for the dependence of the 
diffusion coefficient on friction obtained in detail in the 
Appendix. Section llVI shows our results for the random 
surface potentials. Finally, we conclude with a summary 
and some future directions in Sec. 

II. THE MODEL 

We study the diffusion of a particle in a two- 
dimensional potential V{x, y) of characteristic length 
scale A, in the presence of thermal noise and the associ- 
ated dissipation. The model is embodied in the equations 
of motion, 

d fx y \ 

m * = -fy V [\>\)-rt + S*®> 

where m is the mass of the particle and a dot denotes a 
derivative with respect to t. The generalization to dis- 
tinct length scales X x and X y is straightforward. The 
parameter \x is the coefficient of friction, and the 
are mutually uncorrelated white noises that obey the 
fluctuation-dissipation relation, 



Equations Q can be written in terms of the scaled 
dimcnsionless variables 

where Vq is some measure (e.g. the maximum or the 
mean maximum) of the potential. In terms of these vari- 
ables, and with a dot now denoting a derivative with 
respect to r, 

d 

r x = - j. — V(r x ,r y ) - jr x + ( x (t) 

8 ^ 
r y = —g^-V(r x , r y ) - 7r v + ( v (t), 

where V(r x ,r y ) — V(x/X,y/X)/V is the dimensionless 
potential, and the scaled noise obeys the fluctuation- 
dissipation relation, 

UrMr')) = 2 1 TS ij S(r-r'). (5) 

This scaling serves to stress that there are only two in- 
dependent parameters in this model, the scaled temper- 
ature T and the scaled dissipation 7: 

T = k B T/V 0l 7 = 0AA/roVfc. (6) 

The scaled temperature will be fixed at a value smaller 
than unity so that the thermal fluctuations do not over- 
whelm the potential, but not too much smaller than unity 
so that activated passage over potential barriers is possi- 
ble within a reasonable time. 

The model Q is a "standard" Brownian motion model. 
It relies on well-established ideas of statistical mechanics 
and invokes nothing special about the fluctuations and 
the dissipation. The random forces are normal thermal 
fluctuations, Gaussian and (^-correlated, and the dissi- 
pative forces that accompany the fluctuations are con- 
structed so as to insure thermal equilibration. Neverthe- 
less, we argue that many of the dynamical features of a 
particle evolving under these equations of motion have 
not been investigated until recently 18]. Furthermore, 
we assert that it is not necessary to inject special assump- 
tions such as Levy flights or special memory effects into 
models of surface diffusion, but that these features ap- 
pear naturally from this standard model, thus placing the 
entire range of phenomena observed in surface diffusion 
on an equal common footing. Furthermore, since recent 
experimental and theoretical results have been presented 
for a variety of surfaces, we explore periodic as well as 
random surface potentials, the latter generated accord- 
ing to a given distribution and with a given short-range 
spatial correlation. 

The most straightforward quantity to characterize the 
surface diffusion process is the familiar mean square dis- 
placement of the particle, 



(*0> = 2 l xk B T8 ij 8{t-t'). (2) 



(r 2 (r)) = (Mr) - r^O)] 2 } + ([r y (r) - r,(0)] 2 ). (7) 
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Normal diffusive behavior is characterized by a linear de- 
pendence on time, (r 2 (r)) ~ Dt, where D is the diffusion 
coefficient, while non-diffusive behavior shows a different 
time dependence, (r 2 (r)} ~ t q , with a < 1 for subdiffu- 
sive behavior and a > 1 for superdiffusive motion. Not 
only do we seek to establish the nature of the motion 
as a function of the friction coefficient as measured by 
moments such as the mean square displacement, but we 
seek to establish a theoretical framework that will allow 
us to fully understand and predict the nature of the mo- 
tion under realistic circumstances that may include more 
complex geometries. Toward this purpose we will explore 
the dependence of the diffusion coefficient D on friction 
at large time scales. 

Another characterization of the process at intermedi- 
ate time scales is obtained from the probability density 
function (pdf) P(r, r) of particle displacements r at time 
t. In particular, this pdf reflects the long stretches of 
ballistic motion evident in the intermediate time dynam- 
ics in the low friction regime. Indeed, long ballistic ex- 
cursions lead to a specific behavior of this distribution 
that in turn implies that the velocity of these particles 
remains correlated over considerable time intervals. The 
best way to characterize such correlations is in terms of 
the velocity power spectrum S(cu) — (v(w)v(— u>)). We 
note that this spectrum is directly connected with the 
asymptotic diffusion coefficient of the particle. Accord- 
ing to the Taylor-Kubo formula, D = J* °° C(r)dr, where 
C(t) = (v(t')v(t' + t)) is the velocity-velocity correla- 
tion function. The latter is connected with S(u>) by the 
Wiener-Khinchin relation, S(lu) = J C (r)e %u}r d,T . These 
various observables are studied below for both periodic 
and random surface potentials. 



III. PERIODIC POTENTIAL 
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FIG. 1: A finite portion of the much larger two-dimensional 
periodic potential in which a particle diffuses. 
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FIG. 2: Left: A trajectory for 7 = 1 over t = 20,000 time 
units. Right: A trajectory for 7 = 0.04 over t = 15, 000 time 
units. The period of the potential is A = 4. Note the different 
scales in the two panels. 
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We consider the periodic potential 

^ y )=^ cos(^ + f)co S (^-f), (8) 

which has maxima at positions (nA, mA) and minima at 
((n + |)A, (m + |)A), where n and m are integers. The 
barrier height at the saddle points is Vq. This potential 
is shown in Fig. [I] In Fig. [21 we show examples of trajec- 
tories obtained for two different friction coefficients upon 
numerical simulation of the equations of motion. In our 
simulation we have used A = 4 and m = 1, with a fixed 
scaled temperature T = 0.2 and the scaled dissipation 
as the independent parameter. Note that since the po- 
tential can be generated analytically and the equations 
of motion are continuous, the system is infinite and it 
is not necessary to specify boundary conditions. One of 
the trajectories [panel (a)] is for large friction, and the 
particle is seen to follow what appears to be typical dif- 
fusive motion characterized by short steps of length A 
along the crystallographic directions. The other trajec- 
tory [panel (b)] corresponds to small friction and clearly 



shows the preponderance of long (3> A) tracks along one 
crystallographic direction before turning to another. 

The evolution of the mean square displacement (r 2 ), 
averaged over a set of 5000 particles, is shown in Figs. |3| 
and 0]for several friction coefficients and different initial 
conditions. In Fig. [31 the particles are initially deposited 
in a square of side 2tt around the center of the system ac- 
cording to a Boltzmann equilibrium distribution for the 
positions and for the velocities, while in Fig. 0] they are 
initially located at the center of the system with zero ini- 
tial velocity. Differences in the two cases are observed at 
short times. At very long times the memory of the initial 
condition is lost, and in all cases the motion is diffusive, 
as expected. Regardless of initial condition, for small 
7 and at intermediate times there is clear superdiffusive 
ballistic (a = 2) behavior over several decades of time, 
reflective of the long straight stretches seen in the low-7 
trajectory in Fig. 2] The ballistic behavior is observed 
over a time range of 0(7 _1 ), and one might be tempted 
to conclude that it is therefore a trivial generalization of 
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FIG. 3: Mean square displacement vs time for an ensemble 
of 5000 particles in the periodic potential with a Maxwell- 
Boltzmann initial distribution. 7 = 0.0004 (solid line), 0.004 
(dotted line), 0.04 (dashed line) and 0.4 (dot-dashed line). 



FIG. 4: Mean square displacement vs time for an ensemble of 
5000 particles in the periodic potential with zero velocity, and 
located at the origin. 7 = 0.0004 (solid line), 0.004 (dotted 
line), 0.04 (dashed line) and 0.4 (dot-dashed line). 



the motion of a free damped particle. However, the situa- 
tion here is quite different. The motion is of an ensemble 
of particles whose initial energy distribution is either the 
Boltzmann distribution appropriate to the temperature 
T (Fig. or all of which initially have zero velocity 
(Fig. and the mean square displacement is an average 
over all the particles. The particles with energies lower 
than the barrier height, which includes most of the par- 
ticles, are at first trapped in the original potential well. 
Others (those with a higher initial energy, and those that 
get kicked up sufficiently in energy by the thermal fluctu- 
ations) move out of one well only to be quickly trapped 
again, perhaps into a neighboring well. A few particles, 
those that either start with or acquire the highest veloc- 
ities according to the Maxwell-Boltzmann distribution, 
escape the initial well and move over longer distances 
even as they slow down through friction, until they are 
once again trapped. That this complex ensemble behav- 
ior should give rise to a mean square displacement that 
behaves as shown in Fig. [21 is not self-evident, and in 
fact can only be understood in more detail if we analyze 
the distribution of particle displacements as a function of 
time. We do so later in this section. 

As noted above, for very long times the motion is dif- 
fusive in all cases. The diffusion coefficient D, evaluated 
numerically at long time as 



regime [large 7, see Eq. HA6jl] 



D = lim 



< r 2 (r) > 



4r 



(9) 



is shown in Fig. as a function of the friction 7. The 
solid lines correspond to theoretical predictions evaluated 
in the Appendix. There we show that in the overdamped 



D exp 

7 V T 



(10) 



In the opposite limit [small 7, see Eq. <|A20(l ]. we have 



D 



ttT 

47 



exp 



(11) 



Both predictions fit our numerical results rather well. 
Note that there are no free parameters in these results. 
Note also that while D ~ 7 _1 at both high and low fric- 
tion, the physical reasons are quite different. As seen in 
the Appendix, at high friction the mean square distance 
traveled between trapping events is unity independently 
of 7, and the mean time to escape a well and be trapped 
in another is proportional to 7. At low friction the mean 
square distance is proportional to 7~ 2 and the mean time 
goes as 7 _1 . The ratio of the two is therefore in any case 
proportional to 7" 1 . We note that the low-7 result is not 
applicable to an arbitrary periodic potential |12| . In par- 
ticular, it is correct if, as in our potential, the direction 
in which the particle crosses the saddle point lies along 
the direction of steepest descent. If these do not coin- 
cide, the particle can not simply move along effectively 
one-dimensional channels, and the low-7 behavior of the 
diffusion coefficient may be modified. 

The diffusion coefficient characterizes only the very 
long-time asymptotic dynamics. To characterize the pro- 
cess at intermediate time scales we introduce the proba- 
bility distribution function (pdf) P{r, r) of particle dis- 
placements r at time r. In particular, this pdf reflects 
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FIG. 5: Diffusion coefficient as a function of 7 for the periodic 
potential. 



the long stretches of ballistic motion evident in the inter- 
mediate time dynamics in the low friction regime. The 
pdf is shown in Fig. (top) for 7 = 0.0004 and three 
different time intervals r for an initially equilibrated en- 
semble. For comparison, we also show typical pdf's for 
high damping (7=1, bottom) for the three time in- 
tervals. In the high-7 curves the highest maximum cor- 
responds to no jumps (by far the most likely event at 
short times). The next is associated with jumps to a 
nearest neighbor well, and so on. In contrast, the low- 
7 curves show a very different behavior, with features 
strongly resembling those of a Levy- walk model 0, |2(J : 
a peak at small displacements, a power-law intermediate 
regime, and a side hump at high displacements. Each 
of these is a distinct signature of Levy-walk-like dynam- 
ics, but one must be cautious in the detailed interpre- 
tation of these components, because in our results they 
are not attributable to exactly the same sources as in 
the Levy walk, nor do the exponents fit the Levy walk 
scheme directly. In our system the persistent small dis- 
placement peak is associated with long trapping periods 
during which a particle does not move at all because its 
energy is not sufficient to overcome the barrier. The high 
displacement peak, which moves outward with velocity of 
order unity, is associated with ballistic motion of those 
particles that acquire enough energy to move (and lose it 
very slowly) . Genuine Levy- walk dynamics also exhibit a 
low displacement peak and a superdiffusive peak. In the 
Levy walk, the persistent small displacement peak stems 
from the distribution of excursion lengths, whereas the 
high displacement peak is associated with laminar events 
in which a particle started moving before the observations 
started and are still moving without interruption at time 
t. 

These small- and large-displacement domains are sep- 
arated by a power law behavior. However, here some 



important differences must also be stressed. In the Levy 
walk there is a particular relation between the exponent 
a of the time in the mean square displacement (r 2 ) ~ t a 
and the slope /i of the power law regime of the distri- 
bution, P(r, t) ~ t/t 1 *, namely, a = 4 — fi. Further- 
more, the distribution with the features that we are de- 
scribing is observed in the regime 2 < /z < 3, that is, 
1 < a < 2. In our distribution the intermediate power 
law regime reflects ballistic transport, and yet the slope 
in our power law regime is smaller than unity (approxi- 
mately 0.7). Nevertheless, the qualitative features of our 
distribution track those of the Levy walk. Our side hump 
is strongly broadened whereas the side hump in the Levy- 
walk model is associated with motion at a single constant 
velocity. In our case the velocity varies according to the 
equilibrium Maxwell-Boltzmann distribution. Note that 
the existence of the pronounced side hump moving with 
the velocity of the order of unity reflects the fact that 
the particles performing long steps ("flights") are those 
in the tail of the Maxwcllian velocity distribution. This 
contribution to the pdf in a potential system arises from 
a small subset of particles and is thus to be distinguished 
from that of a typical underdamped free Brownian par- 
ticle plj . In summary, in our underdamped system the 
slower particles are trapped and localized around a po- 
tential minimum. They do not contribute to transport 
and lead to the first peak in the pdf P(r, r). Parti- 
cles with higher energies begin moving essentially bal- 
listically, thus contributing to the side hump of the pdf 
that moves toward the right with increasing times in the 
first panel in Fig. These particles eventually become 
trapped due to friction, and these trapping events lead to 
the power law portion of the pdf. Fewer and fewer parti- 
cles (only those in the ever smaller highest-velocity tails 
of the distribution) continue moving ballistically. The 
side hump thus becomes narrower and lower with time. 
At long times, r 3> I/7, both the long trapping and bal- 
listic features are of course no longer present as diffusive 
motion dominates the behavior. 

We emphasize that the behavior described above is uni- 
versal and is not strongly dependent on whether the ini- 
tial conditions are equilibrium or zero-velocity ones (al- 
though the exact heights of the peaks do). In Fig. [7| the 
velocity probability distributions of our system at finite 
times are compared to a Maxwell distribution. When 
we implement Maxwell-Boltzmann initial conditions, the 
overall velocity distribution obtained in our simulations 
by time-sampling is purely Maxwellian at all times. For 
the noncquilibrium initial conditions where we start with 
zero velocity at the saddle point of the potential, some 
nonequilibrium transient phenomena are observed in the 
small deviations from the Maxwell form. We also note 
that the velocity distribution thermalizes at times of 
the order of 107 -1 . The times in Fig. [7| correspond to 
t = 27- 1 and 47- 1 . 

The long, ballistic excursions that lead to the power 
law and side hump contributions to the particle dis- 
placement distribution function imply that the veloci- 
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FIG. 6: Log-log plot of P(r, r) for 7 = 0.0004 (top) and for 
7 = 1 (bottom) and three different values of time intervals : 
t — 20 (squares), 100 (open circles), 1000 (triangles). Note 
the strong difference in scales. 
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FIG. 8: S(u) for 7 = 0.004, 0.04, 0.4 and 4 in decreasing 
order. 



ties of these particles remain correlated over consider- 
able time intervals. This aspect is most directly char- 
acterized by a study of the velocity power spectrum 
S(uj) = (v(lj)v(—u))). In the case of diffusive motion 
one can extract the diffusion coefficients as, 

D = 1 (S xx (lj = 0) + Syy(u = 0)) , (12) 

where Su(u>) is the power spectrum of the cartesian com- 
ponent i of the velocity vector. The behavior of S(w) for 
different values of the friction coefficient 7 is shown in 
Fig. |SJ The structure of the power spectrum mirrors the 
underlying dynamics. Thus, an evident feature of S(oj) 
is the existence of pronounced peaks at the frequency 
ujq = ny/2 of small oscillations in one well performed by 
particles trapped in a well. At low frequencies, wo < 1, 
a power-law growth of S(u>) is observed. This is reflec- 
tive of the persistent time correlations associated with 
ballistic excursions. The slope ~ 1 of this growth corre- 
sponds to C(t) ps const and to the ballistic growth of the 
mean square displacement. At even smaller frequencies 
the power spectrum crosses over to S(co) = const = 4£>, 
indicating full decorrelation and emergence of pure diffu- 
sion. 



FIG. 7: Probability distribution of the velocity for two differ- 
ent initial conditions compared with a Maxwell distribution 
(solid curve) for the same temperature. Triangles: distribu- 
tion at t — 50 for zero velocity initial condition. Circles: 
distributions at r = 100 for zero velocity (solid circles) and 
Maxwell (open circles) initial distributions (7 = 0.04). 



IV. RANDOM POTENTIAL 

Surfaces arc usually not completely crystalline or reg- 
ular because of the presence of vacancies, defects, and 
other types of disorder. One type of surface disorder is 
represented by a random potential with spatial correla- 
tions that can model the presence of some spatial finite- 
ranged order. An algorithm that can be used to generate 
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FIG. 10: Left: A trajectory for 7 = 0.1 over t = 250000 
time units. Right: A trajectory for 7 = 0.0001 over t — 2000 
time units. The characteristic length scale of the potential is 
A' = 4. Note the different scales in the two panels. 



FIG. 9: Random potential with the same average potential 
height as in the periodic case. 

surfaces with any given spatial correlation has been pre- 
sented in Refs. |22l 1231 124^. We have implemented this 
algorithm for a random potential surface with a Gaus- 
sian distribution and exponential correlation function. In 
terms of the vectors x = (x, y) and x' = (x' , y'), our po- 
tential surface has the correlation property, 

(nx)l/(x'))= S (x-x'), (13) 

with 

3 (x-x') = ^ e -'— 'I 2 /-' 2 , (14) 

which is parameterized by the intensity e and the char- 
acteristic length A' that we take to be the same as the 
length scale in the periodic potential (A' = 4). Note that 
the average potential height is a combination of the pa- 
rameters. A typical surface generated with this algorithm 
and of average potential height Vq equal to that of the 
periodic potential is shown in Fig. El In our particular 
case, we have done our simulations on a square grid of 
size L = 4096 with Sx — 1.0. Even though our potential 
is generated on a grid, extrapolation of forces can be used 
to calculate the force at any point, which is a necessary 
step in the dynamical simulation. We have used the pa- 
rameter e — 100. In this case we must generate finite 
systems. Rather than impose boundary conditions, we 
stop each simulation sequence whenever the first particle 
reaches a boundary. 

Typical trajectories associated with this potential at 
high and low friction are shown in Fig. EH Again, the tra- 
jectories are entirely different, that associated with high 
friction consisting of extremely short steps followed by 
random changes in direction while that associated with 
low friction again shows a pronounced directional persis- 
tence. 

An analysis of the exponents in the mean square dis- 
placement formula shows the entire range of behaviors 
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FIG. 11: Mean square displacement for a particle in the ran- 
dom potential with Maxwellian initial velocity distribution, 
for 7 = 0.0001 (solid line), 0.001 (dotted line), 0.003 (dashed 
line) and 0.008 (dot-dashed line) 



from subdiffusive to superdiffusive with changing fric- 
tion. In Fig^] we show (r 2 (r)), averaged over 5000 
particles, as a function of time for several values of 7. 
A detailed analysis of these trajectories presents a rich 
phenomenology of possible different regimes. In the over- 
damped regime we observe clearly subdiffusive behavior 
(a < 1), already known from overdamped simulations 
presented earlier j24[. An interesting outcome is the su- 
perdiffusive (a > 1) behavior seen for small values of 7. 
This again is a signature of the Levy-walk-type behavior. 
The exponent a as a function of 7 is plotted in Fig. 1121 
Each exponent is obtained as an average over 500 par- 
ticles in each of 20 different realizations of the random 
potential, the error bars indicating the standard devia- 
tions over these statistics. For large 7 our simulations run 
for a long time and the value of a has been calculated on 
the basis of the behavior of (r 2 ) over the last time decade 
of the simulation for each run. In the medium and low 
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FIG. 12: Exponents a versus friction coefficient 7. Each 
point is obtained from an average over 20 realizations of the 
random potential and 500 particles per realization. The error 
bars indicate the standard deviation over these statistics. 



damping cases the simulations do not extend over such 
long times (since each simulation ends when a particle 
reaches the system boundary). The standard deviations 
are larger for small 7. In part this is due to the fact that 
for small 7 it is more likely for a particle to reach the 
boundary of the system during the time of a run. Not 
only are the run times shorter for smaller 7, but they 
are also more broadly distributed as 7 decreases. In any 
case, the trend in the behavior of a as a function of the 
friction is clear. An unexpected result is the absence of a 
diffusive regime (a « 1) except in a small range of values. 

The behavior of S(u>) for the random potential is seen 
in Fig. E| for several values of 7. For small 7 the spec- 
tral density grows with decreasing frequency, indicating 
superdiffusion, while at large 7 it decays with decreasing 
to, which is a signature of subdiffusion. For the interme- 
diate value 7 = 0.003 (close to the value at which a « 1 
is observed in Fig. I12JI it flattens, as appropriate for dif- 
fusive behavior at long times. 

In Fig. ^] the pdf P(r, r) for the random potential is 
plotted. This pdf looks very different from the one for 
a periodic potential. Nevertheless, for the case of small 
friction it still exhibits some features of Levi-walk-likc 
dynamics. The main difference here is the absence of the 
pronounced central peak due to trapping in the nearest 
potential well. This trapping is still evident in the high 
friction trajectory in Fig. 1101 However, for small friction 
the corresponding peak in the pdf is smoothed out by 
the broad distribution of depths of the sampled minima. 
Instead, the distribution resembling a smoothed form of 
the genuine ballistic Levy- walk pdf appears, character- 
ized by the central power law dip and side peak. The 
side peak moves outward with a constant velocity and 
becomes broadened by scattering and trapping. At long 
times the pdf develops the tent-like central peak typical 



FIG. 13: S(u>) for 7 = 0.0001, 0.003 and 0.08 (decreasing or- 
der) in the random potential with Maxwellian initial velocity 
distribution. 




FIG. 14: Log-log plot of P(r,r) for 7 = 0.0001 (top) and 
for 7 = 0.08 (bottom) and various different values of time 
intervals: r = 10 (squares), 20 (open circles), 100 (triangles), 
1000 (rhombuses). 



of trapping. 
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V. CONCLUSIONS 

We have investigated the behavior of an ensemble of 
particles in a two-dimensional potential subject to ther- 
mal fluctuations described by ordinary Langevin dynam- 
ics. Our purpose has been to determine whether such 
an ordinary description is sufficient to produce the entire 
range of behaviors, from subdiffusive to superdiffusive, 
that has been observed in the motion of molecules on 
surfaces [J|lH&!EHSS[ll[H0. Indeed, we 
found that it is not necessary to add any further assump- 
tions to the model to observe the full range of behaviors, 
at least on intermediate time scales. For a given surface 
potential one can introduce scaled parameters such that 
the only control parameters in a periodic potential are 
the temperature (which we have held fixed in this anal- 
ysis) and the dissipation parameter 7. In the case of a 
random potential we introduced an additional parameter 
(which we also hold fixed in this analysis) to character- 
ize the intensity of the spatial correlation function of the 
potential variations. 

For the periodic surface potential we found that in the 
underdamped regime the motion of the particles includes 
a ballistic range that can extend over many decades of 
time. The probability distribution function of the dis- 
placement of the particles exhibits a structure reminis- 
cent of that associated with Levy walks. The pdf has a 
peak at short distances that arises from those particles 
that are trapped in their original well. The pdf also ex- 
hibits a hump that moves outward linearly with increas- 
ing time, and that is due to those particles whose energy 
is sufficiently high for ballistic motion over the time of ob- 
servation. As these particles in turn get trapped, the am- 
plitude and width of this hump decrease. Particles that 
start out moving ballistically eventually lose enough en- 
ergy through damping to get trapped, and these progres- 
sive trapping events give rise to a power law intermediate 
regime connecting the trapped particle peak and the bal- 
listic particle hump. These behaviors are not observed in 
the case of high dissipation. The long stretches of essen- 
tially ballistic motion of the more energetic particles also 
make themselves apparent in the velocity power spec- 
trum. We have ascertained that the Levy-like features 
in the low damping case arise from a subset of energetic 
particles by confirming that throughout the evolution the 
velocity distribution of particles remains Maxwellian if 
the initial distribution is Maxwellian. If initially the par- 
ticles have zero velocity, the distribution rather quickly 
settles into an essentially Maxwellian distribution. Thus, 
we stress again that in our model these features arise not 
because the average particle is subjected to any esoteric 
fluctuations, but rather from the fact that the energetic 
particles in an ordinary Maxwell distribution can move 
over long distances when the system is underdamped. At 
long times (low frequencies) the motion of the particles 
eventually becomes diffusive. This is most clearly seen 
in the behavior of the mean square displacement (which 
eventually grows linearly with time), and in the plateau 



of the velocity power spectrum at low frequencies. 

While most of the results reported above are obtained 
numerically, we are able to calculate the diffusion coef- 
ficient D that describes the long-time behavior of the 
system analytically for high damping (a well-known re- 
sult) and also for low damping. The agreement of our 
analytic results for D with numerical ones is good over 
the entire range of damping coefficients. 

The case of random potential gives rise to a large va- 
riety of regimes ranging from superdiffusion to subdif- 
fusion. While the superdiffusive behavior is probably a 
transient, just as in the periodic potential, subdiffusion 
due to multiple trapping is known to be a true asymp- 
totic behavior when the temperature is low. However, in 
the case of a random potential the restrictions imposed 
by our numerical procedure are tighter, so that the true 
asymptotic behavior is difficult to reach. At short times 
and for weak damping the Levy-walk-like features aris- 
ing from long stretches of motion of energetic particles 
are still evident in our simulations. For strong damping 
this regime is pertinent only for the shortest times and 
then crosses over to that typical of subdiffusive, disper- 
sive transport. 
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APPENDIX A: DIFFUSION COEFFICIENT FOR 
HIGH AND LOW FRICTION IN A PERIODIC 
POTENTIAL 

We start with the assumption that a particle that has 
sufficient energy to move away ( "escape" ) from a poten- 
tial well will in general preferentially move along direc- 
tions of lowest potential barriers, that is, along the r x 
direction at r y — 1/2 (or an odd multiple thereof) or 
along the r y direction with r x — 1/2 (or an odd multiple 
thereof). Choosing the former, we then have the particle 
moving along a line in the periodic potential 

V(r x ) = V(r x ,r y = 1/2) = -sin 2 ^,. (Al) 

Only those particles with energy E > can move. As 
they move they lose energy through dissipation, until 
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they become trapped. We calculate the diffusion coef- 
ficient for the ensemble of particles according to the for- 
mula 



D=(l 2 )/(2r). 



(A2) 



Here (I 2 ) is the mean square displacement from the ini- 
tial well to the well associated with the next trapping 
event, r is the mean time for such a journey, and the 
numerical factor in the denominator is 2 instead of 4 [cf. 
Eq. |JSJ|] because the motion along each stretch is one- 
dimensional. The average is over all the particles in the 
ensemble, which includes those that initially have suffi- 
cient energy to move away from a potential well and also 
those that do not, as dictated by the Maxwell-Boltzmann 
distribution. For low-energy particles to be able to jump 
out they must acquire enough energy from the thermal 
fluctuations. 

The mean time r in general includes the time for a 
particle to acquire enough energy (E > 0) to move out 
of a well, and the time to become trapped again in a 
new well once it has escaped. We can calculate r using 
arguments for the calculation of transition rates of the 
Kramers problem 25] by associating it with the inverse 
of a transition rate that consists of two contributions: 



1 

2k 



2fc* 



(A3) 



One contribution is the transition state theory rate, 



k tst = £0 e -i/r 



(A4) 



where loq is the angular frequency at the well bottom, 
which in our dimensionless variables and for our potential 
is given by uj = [V"(l/2)] 1 / 2 = tt\/2. The other is the 
transmission factor re, which depends in a non-trivial way 
on the damping. 

Consider first the case of high friction. In the over- 
damped regime (large 7) we assume that (I 2 ) = 1 , an d 
the transmission factor re has the well-known form 25], 



1 



1/2 



(A5) 



where lo 2 = V"(r&) = 2n 2 and r b is the position of the 
maximum of the potential. Thus the diffusion coefficient 
in this regime is 



D = 



1 

27 



2n' 



1/2 



-l/T 
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(A6) 



Now consider the diffusion coefficient in the small-7 
limit or infradamped regime. In this regime a more elab- 
orate analysis is required to calculate (I 2 ). The equations 
of motion of a particle in the potential Eq. (|A1|) can be 



expressed in the standard momentum-position form as 
the pair, 



r x = V 

p = -7p - V'(r x ) + C(r). 



(A7) 



In the small-7 limit the energy variation of a particle is 
slow, and it is more convenient to rewrite the dynamical 
equations in terms of the displacement and the energy 
E=p 2 /2 + V(r x ): 



E 



{2[E-V(r x )}} 1/2 

-2 7 [E-V{r x )]+{2[E- 



(A8) 



V(r x )}} 1/2 m. 



(A9) 



Two quantities needed for our estimate of the diffusion 
coefficient can be calculated from these dynamical equa- 
tions. One is the time t(E) that it takes a particle of 
energy E to traverse a unit distance (the scaled spatial 
period). The other is the energy AE lost by the par- 
ticle when it traverses this distance [2(|. To calculate 
t(E) we assume that the energy E remains fixed during 
the traversal, so that we can simply integrate Eq. (|A8|) 
for constant E. With the potential IjAljl one obtains a 
standard integral: 



-(E) = f 
Jo 



dr x 



V2X 



-.K 



1 



10 y/2 [E-V{r x )\ irVT+E \ V 1 + E ) 

(A10) 

where K, is the complete elliptic integral of the first 
kind [2^]. For E <C 1 the elliptic integral can be ap- 
proximated by a logarithmic leading term, which leads 
to, 



r(E) 



tta/2 E 



(All) 



Next, we estimate the energy loss that occurs during 
this traversal. This has been calculated in many ways, 
but a particularly transparent argument [26j is to ne- 
glect the fluctuations in Eq. i|A9|l (since particles above 
the barrier primarily lose energy) and integrate over a 
time interval t(E). The r x dependence in Eq. (|A9|) even 
after this approximation still poses a problem because it 
leads to incomplete elliptic integrals amenable only to nu- 
merical integration. Since r x changes rapidly compared 
to E, it is reasonable to perform an average of the r x - 
dependent term over a traversal from to 1. This leads 
to the approximate equation, 



E = -2 1 F{E) 



(A12) 



where, 



with, 



F(E)= f dr x P E {r x )[E-V{r x )] (A13) 
Jo 



P E {r x ) 



1 



t{E) [E-V(r x )Y 



(A14) 
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The result is, 




(A15) 



where £ is the complete elliptic integral of the second 
kind |2jj • Retention of terms to leading order in E and 
AE/E leads to 



AE : 



lV8 



(A16) 



which is independent of E. 

From the result for AE we can now calculate the mean 
square distance (I 2 ) traveled by a particle that starts with 
an energy above the barrier. Since the particle loses en- 
ergy AE in each traversal of a unit distance, the distance 
traveled by such a particle of initial energy E before being 
trapped is, 



v ' AE 



nE 



(A17) 



Averaging over a Boltzmann distribution of energies, we 
obtain, 



(I 2 ) = (k B T)- 



dE l 2 (E)e 



-E/T _ 



4 7 2 



(A18) 



When a particle of positive energy E moves from one 
potential well to the next, it does so in a time t(E) and 
it loses energy AE in the process. If its energy is now 
below the barrier, it will be trapped. If not, it will keep 
moving, and will take a time t(E—AE) to get to the next 
potential well. The particle moves on until its energy is 
too low to continue moving. The energies have a Boltz- 
mann distribution, so different particles will take differ- 
ent times to be trapped both because they are moving 
progressively more slowly and because they get trapped 
at different times in their journey. Taking these effects 
into account, one finally arrives at the following expres- 
sion for the transmission factor k 1261: 



tanh 



\AE\ 
2T 



—AE 
2T ' 



(A19) 



where the last expression holds for very low damping. We 
thus finally obtain the following result for the diffusion 
coefficient in this regime: 



D 



ttT 

47 



-l/T 



(A20) 
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